open Glob

let create = Array.make_matrix

let init nx ny f =
    Array.init nx (fun ix -> Array.init ny (f ix))

let make_zero m = 
    let nx = Array.length m and ny = Array.length m.(0) in
    for x = 0 to nx-1 do 
        for y = 0 to ny-1 do 
            m.(x).(y) <- 0.
        done
    done

let print ?(diagonal=false) arr =
  let l = Array.length arr in
    for i = 0 to l - 1 do
      (* print only half and diagonal of this symmetric matrix *)
      print10 (Printf.sprintf "%i |" (i + 1)) 
        (if diagonal
         then (Array.sub arr.(i) 0 (i + 1))
         else arr.(i))
    done

let map f m =
    Array.map (fun a -> Array.map f a) m

let map_ij f m =
    Array.mapi (fun i a -> Array.mapi (fun j elt -> f (i,j) elt) a) m

let iter f m =
    Array.iter (fun a -> Array.iter f a) m

let iteri f m =
    Array.iteri (fun i a -> Array.iteri (f i) a) m

let iterk f m =
    Array.iteri (fun i a -> Array.iteri (fun j v -> f (i,j) v) a) m

let get m (i,j) = m.(i).(j)

let set m (i,j) v = m.(i).(j) <- v

let transpose a =
    let nx = Array.length a and ny = Array.length a.(0) in
    let a'T = Array.make_matrix nx ny 0. in
    for kx = 0 to nx - 1 do
        for ky = 0 to ny - 1 do
            a'T.(ky).(kx) <- a.(kx).(ky)
        done
    done;
    a'T

(*
let sum arr =
     Array.fold_left (fun sum array ->
         Array.fold_left (
             fun acc c -> if not_nan c then acc+.c else acc
         ) sum array
     ) 0. arr
let minmax arr=
    let min,max=ref arr.(0).(0), ref arr.(0).(0) in
    for kx=0 to Array.length arr - 1 do
        for ky=0 to Array.length arr.(kx) - 1 do
            let a=arr.(kx).(ky) in
            if a > !max then max := a;
            if a < !min then min := a
        done
    done;
    !min, !max
*)

let fold_left f init matrix =
    Array.fold_left 
      (fun acc arr -> Array.fold_left (fun acc a -> (f acc a)) acc arr) 
      init matrix

include Glob.Make_safe_float_op(struct
    type t = float array array type place_info = int*int
    let fold = fold_left let iterk = iterk
end)

let fold_lefti f init matrix =
    let acc = ref init in
    Array.iteri (fun i arr ->
        Array.iteri (fun j elt -> 
            acc := f !acc (i,j) elt
        ) arr
    ) matrix;
    !acc 

